#在服务器上运行代码
setwd("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20210324")
library(openxlsx)
filea=read.csv("all.FDR.sig.at.least.one.add.direction.same.diff.csv",head=T)
filea$id=paste(filea$Chr,filea$Start,sep = ":")
filea1=filea[filea$FDR.sig>1,]
file=read.table("53K.add.GWAS.eQTL.DEG.motif.for.analysis.txt",header=T,sep="\t")
file=file[!duplicated(file$id),]	#对TF去重
file1=file[file$pattern.not.rm.dupl.num.DC>1,]
file2=file1[file1$BF_in_DC>1,]
file3=file1[file1$BF_in_DC>10,]

asm2=read.table("220520ASMs_anno.hg19_multianno.csv",head=T,sep=",")	#220K
asm2$unitID=paste(asm2$Chr,asm2$Start,asm2$Ref,asm2$Alt,sep=":")
asm=read.csv("869727.all.snp.vaf.up.down.20210321.csv",head=T)
asm=asm[asm$unitID %in% asm2$unitID,]
asm=asm[as.character(asm$unitID) %in% as.character(filea$unitID),]
#asmdata=filea1[filea$unitID %in% asm$unitID,]
asm1=tidyr::separate(asm,unitID,into=c("chr","pos","ref","alt"),sep=":")
asm1$id=paste(asm1$chr,asm1$pos,sep=":")

###control
con=read.table("all.heter.snp",head=T,sep="\t")
con=con[!con$id %in% as.character(filea$id),]		###remove the ASH
conid=as.character(con$id)
conid.random=sample(conid,150000,replace = FALSE)



#LIBD eQTL
anno=fread("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207/BrainSeqPhaseII_snp_annotation.txt",head=T,sep="\t")
anno=data.frame(snp=anno$snp,chr_hg19=anno$chr_hg19,pos_hg19=anno$pos_hg19)
anno$id=paste(anno$chr_hg19,anno$pos_hg19,sep=":")

filea=merge(filea,anno,by="id",all.x=T)
filea1=merge(filea1,anno,by="id",all.x=T)
file=merge(file,anno,by="id",all.x=T)
file1=merge(file1,anno,by="id",all.x=T)
file2=merge(file2,anno,by="id",all.x=T)
file3=merge(file3,anno,by="id",all.x=T)
asm1=merge(asm1,anno,by="id",all.x=T)

library(data.table)
sel=list.files(pattern="BrainSeqPhaseII",path="/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207",full.names=T)[-13]
tiss =c("DLPFC","HIPPO","Interaction")
regions =c("exon","gene","jxn","tx")

conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
con1=merge(con1,anno,by="id",all.x=T)

						####################研究ASH在各个组织各染色体位置（"exon","gene","jxn","tx"）的调控作用，这是比较细分的。

for(randomid in 1:3){
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
con1=merge(con1,anno,by="id",all.x=T)

result=data.frame(matrix(,nrow = 1,ncol =6))
names(result)=c("Term","Tissue","regions","overlap.num","OR","P.value")
result=result[-1,]

for(i in tiss){
for(j in regions){
fn=paste("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207/BrainSeqPhaseII_eQTL_FDR1perc_",i,"_",j,".txt",sep="")
ftmp=fread(fn,header=T,sep="\t")

c=length(intersect(as.character(con1$snp),as.character(ftmp$snp)))
d=150000-c
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
rt[1,]=c("control",i,j,c,NA,NA)
result=rbind(result,rt)

a=length(intersect(as.character(filea$snp),as.character(ftmp$snp)))
b=dim(filea)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("117012",i,j,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(filea1$snp),as.character(ftmp$snp)))
b=dim(filea1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("61725",i,j,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file$snp),as.character(ftmp$snp)))
b=dim(file)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("53425",i,j,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file1$snp),as.character(ftmp$snp)))
b=dim(file1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("8544",i,j,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file2$snp),as.character(ftmp$snp)))
b=dim(file2)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("807",i,j,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file3$snp),as.character(ftmp$snp)))
b=dim(file3)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("200",i,j,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(asm1$snp),as.character(ftmp$snp)))
b=dim(asm1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =6))
names(rt)=c("Term","Tissue","regions","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("13649.117KASH.22OKASM",i,j,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)
}

}
result$id=paste(result$Tissue,result$regions,sep="_")
rtn =paste("tissue.region.LIBD.eQTL.enrich.random",randomid,".csv",sep="")
write.csv(result,rtn,quote=F,row.names=F)
}

						####################研究ASH在2个脑组织（包含3种方法："DLPFC","HIPPO","Interaction"）的调控作用


for(randomid in 1:3){
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
con1=merge(con1,anno,by="id",all.x=T)

for(i in tiss){

result=data.frame(matrix(,nrow = 1,ncol =5))
names(result)=c("Term","Tissue","overlap.num","OR","P.value")
result=result[-1,]

ft=data.frame(matrix(,nrow = 1,ncol =3))
names(ft)=c("snp","Symbol","Type")
ft = ft[-1,]

for(j in regions){
fn=paste("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207/BrainSeqPhaseII_eQTL_FDR1perc_",i,"_",j,".txt",sep="")
ftmp=fread(fn,header=T,sep="\t")
ftmp=data.frame(snp=ftmp$snp,Symbol=ftmp$Symbol,Type=ftmp$Type)
ft =rbind(ft,ftmp)
}
c=length(intersect(as.character(con1$snp),as.character(ft$snp)))
d=150000-c
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
rt[1,]=c("control",i,c,NA,NA)
result=rbind(result,rt)

a=length(intersect(as.character(filea$snp),as.character(ft$snp)))
b=dim(filea)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("117012",i,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(filea1$snp),as.character(ft$snp)))
b=dim(filea1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("61725",i,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file$snp),as.character(ft$snp)))
b=dim(file)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("53425",i,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file1$snp),as.character(ft$snp)))
b=dim(file1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("8544",i,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file2$snp),as.character(ft$snp)))
b=dim(file2)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("807",i,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file3$snp),as.character(ft$snp)))
b=dim(file3)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("200",i,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(asm1$snp),as.character(ft$snp)))
b=dim(asm1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =5))
names(rt)=c("Term","Tissue","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("13649.117KASH.22OKASM",i,a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

rtn =paste("tissue.",i,".LIBD.eQTL.enrich.random",randomid,".csv",sep="")
write.csv(result,rtn,quote=F,row.names=F)
}
}


						####################研究ASH在总的LIBD eQTL调控作用


for(randomid in 1:3){
conid.random=sample(conid,150000,replace = FALSE)
con1=con[con$id %in% conid.random,]
con1=merge(con1,anno,by="id",all.x=T)

result=data.frame(matrix(,nrow = 1,ncol =4))
names(result)=c("Term","overlap.num","OR","P.value")
result =result[-1,]

ft=data.frame(matrix(,nrow = 1,ncol =3))
names(ft)=c("snp","Symbol","Type")
ft = ft[-1,]
for(i in tiss){

for(j in regions){
fn=paste("/mnt/md1200/6/yjp/5hmc_analysis_hg19_new/20201207/BrainSeqPhaseII_eQTL_FDR1perc_",i,"_",j,".txt",sep="")
ftmp=fread(fn,header=T,sep="\t")
ftmp=data.frame(snp=ftmp$snp,Symbol=ftmp$Symbol,Type=ftmp$Type)
ft =rbind(ft,ftmp)
}
}
c=length(intersect(as.character(con1$snp),as.character(ft$snp)))
d=150000-c
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
rt[1,]=c("control",c,NA,NA)
result=rbind(result,rt)

a=length(intersect(as.character(filea$snp),as.character(ft$snp)))
b=dim(filea)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("117012",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(filea1$snp),as.character(ft$snp)))
b=dim(filea1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("61725",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file$snp),as.character(ft$snp)))
b=dim(file)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("53425",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file1$snp),as.character(ft$snp)))
b=dim(file1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("8544",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file2$snp),as.character(ft$snp)))
b=dim(file2)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("807",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(file3$snp),as.character(ft$snp)))
b=dim(file3)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("200",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

a=length(intersect(as.character(asm1$snp),as.character(ft$snp)))
b=dim(asm1)[1]-a
rt=data.frame(matrix(,nrow = 1,ncol =4))
names(rt)=c("Term","overlap.num","OR","P.value")
tmp=fisher.test(matrix(c(a,c,b,d),nrow = 2,byrow = T))
rt[1,]=c("13649.117KASH.22OKASM",a,tmp$estimate,tmp$p.value)
result=rbind(result,rt)

rtn =paste("LIBD.eQTL.enrich.random",randomid,".csv",sep="")
write.csv(result,rtn,quote=F,row.names=F)
}





